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ABSTRACT 

The observed cooling rate of hot gas in clusters is much lower than that inferred from 
the gas density profiles. This suggests that the gas is being heated by some source. 
We use an adaptive-mesh refinement code (Flash) to simulate the effect of multiple, 
randomly positioned, injections of thermal energy within 50 kpc of the centre of an 
initially isothermal cluster with mass M200 = 3 x 10^'' Mq and kT — 3.1 keV. We 
have performed eight simulations with spherical bubbles of energy generated every 
10* years, over a total of 1.5 Gyr. Each bubble is created by injecting thermal energy 
steadily for 10^ years; the total energy of each bubble ranges from 0.1-3x10^" erg, 
depending on the simulation. We find that 2x10^° erg per bubble (corresponding to a 
average power of 6.3 x 10^^ erg s~^) effectively balances energy loss in the cluster and 
prevents the accumulation of gas below kT = 1 keV from exceeding the observational 
limits. This injection rate is comparable to the radiated luminosity of the cluster, and 
the required energy and periodic timescale of events are consistent with observations 
of bubbles produced by central AGN in clusters. The effectiveness of this process 
depends primarily on the total amount of injected energy and the initial location of 
the bubbles, but is relatively insensitive to the exact duty cycle of events. 
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1 INTRODUCTION 

Gas cooling at the centre of a cluster halo is an inherently 
unstable process: cooling increases the gas density, which 
in turn enhances the cooling rate. X-ray observations show 
that the cooling time of g as in most cluster core s is less 
than the Hubble time (e.g. iFabian fc Nulsenlll977^ . Unless 
cooling is balanced by some form of heating, gas will flow 
into the cluster cent re at rates up to ~ 1000 MQyr~^ (e.g. 
iPeterson era!] 1200 j) . What happens to the cooled gas is 
unclear; it could be consumed by star formation, or lead 
to a reservoir of low temperature (< 1 keV) material in 
the core. Although this may be the fate of a fraction of 
the gas, there are indications that most gas in fact does 
not follow either route. First, the star formation rates in 
the central galaxies of clus ters are much lower than the 
inferred mass inflow rates dO'Connell fc McNamaral Il989l : 
Ijohnstone. Fabian, fc Nulsenlll987ll . rarely approaching ~ 
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100 M0yr~^. To estimate the star formation rate in a typi- 
cal cluster, we can average o ver th e large sample of clusters 
studied bv ICrawford et all (Il999ll . This suggests that the 
typical star formation rate is less than 10 MQyr~^. Secondly, 
we can compare the mass deposition rate s with observ ations 
of the molecular gas content of clusters. lEdgel Jiooi finds 
molecular gas masses ranging from lO'' to 2 X lO" Mq, with 
an average of 2.6 X 10^" Mr Ti. Assuming a gas consumption 
timescale of 10 years fsee lEdgeir2001ll . this implies a de- 
position rate that may be as high as 200 MQyr~^ in a few 
clusters, bu t is ^ 30 Mr^yr~^ on aver age. Similar limits are 
obtained bv lSalome fc Combed i2003l) . Finally, recent spec- 
troscopic X-ray observat ions show no evidence for significant 
gas cooling bel ow 1 keV ([Kaastra et all200ll : lPeterson et alJ 
[2001; Tamura erai] 1200 ill , and observations of molecular 
and neutral material reveal that the amount of cold gas in 
clusters of galaxies is also much less than expected from the 
in tegrated cooling flow rate (typically less than 30 M^yr~^ 
- lEdge et al.ll2002l : lidge fc Frayeil20oilSalome fc Combed 
l2003h . 
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This cooling-flow paradox has led many authors to 
investigate mechanisms to quench gas cooling. Observa- 
tions of merging clusters show little evidence for cool- 
ing flows, which suggests that the merger process might 
be implicated in disrupting cooling. Recent simulations 
have sho wn that sub-halo merging can indeed heat up gas 
iBurns e t al. 2003); however, the amount of cold gas pro- 
duced in these simulations is still too large compared to 
that observed, leading to the conclusion that additional 
heating processes must be involved. Several alternatives 
have been proposed, including energy injection from ra- 
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viscous dissipation of sound waves ([Fabim^^l 12003. , 
lRuszkowski^ _^ruggenj_fc" 3egelnian||200g) and the rmal con- 
duction JVoigt fc FabianI 120031: iDolag et al.ll2004r) . Each of 
these has advantages and disadvantages. In particular, it is 
difficult to balance cooling, with its density dependence, 
with heating processes that typically scale as p. Since cooling 
and heating can then balance only at one particular density, 
some of these feedback mechanisms may lead to unstable 
solutions with some regions of the cluster being efficiently 
heated while others continue to cool catastrophically. 

AGN are particularly promising candidates for balanc- 
ing cooling, given their potentially large energy reservoir 
jTabor fc Binnevlll993l: iBinnev fc Tabojll995l: iBower et alJ 



20011) . In most numerical simulations of this process to 
date, energy injection produces bubbles at high tem- 
perature and low density that, after a short expansion 
phase , gain momentum by buoyanc y (|BruggerL.fc.JCaiser 
20011: iQuifis. Bower, fc Baloghl 1200 it iBasson fc Alexander 
2003r . This mimics the effect of a jet which rapidly loses 
its co Uimation, as often observed in local clusters jEilekl 
|20o3). In other simulations, gas is injected at high veloc- 
ity to mimic jets which retain their large-scale cohere nce 
jRevnolds. Heinz, fc Begelm an"2002': tomma et alJl2003t) . In 
this paper, we consider the first mechanism and show that 
such bursts of localised energy can induce convection of the 
intra-cluster medium (ICM), which leads to a quasi-stable 
cluster configuration and reduces the average mass deposi- 
tion rate to within observational limits. We discuss whether 
the required energy and duty cycle of AGN activity are com- 
patible with observational limits. This paper is organised as 
follows. In Section|5|we describe the setup of our simulations 
and discuss the parameters explored. Results are presented 
in Section 13 and, in Section |31 we compare the required en- 
ergy and duty cycle with observational constraints. Finally, 
our conclusions are summarised in Section |S| 



2 SIMULATIONS 

In order to concentrate our numerical resolution on the in- 
teractions between the cooling material and the buoyant hot 
bubbles, our simulations use a fixed external gravitational 
potential and neglect the self-gravity of the gas. This al- 
lows us to run the simulation over cosmologically significant 
timescales. We model the input of energy as a cycle of en- 
ergy injection and quiescent phases. In the active phase a 
small region of the gas is heated so that it expands to form 
a bubble that rises out of the confining potential. 



name 


E 

lO^o erg 


Injected/Radiated Power 

Ei < Er >ioo 
10*4 erg s-i lO''-* erg s"! 


A 





U 


no cooling 


SO.O 








28 


SO.l 


0.1 


0.32 


18 


S0.3 


0.3 


0.95 


18 


S0.6 


0.6 


1.9 


13 


Sl.O 


1.0 


3.2 


13 


SI. 5 


1.5 


4.7 


8.7 


S2.0 


2.0 


6.3 


7.9 


S3.0 


3.0 


9.5 


6.8 



Table 1. Summary of the nine simulations performed. The dif- 
ferent simulations are referred to as Sn, where the energy i? of a 
single bubble is n X 10^^ erg. Ei is the mean energy injection rate 
over a duty cycle, and < Er >ioo is the mean emitted energy 
rate, averaged over the last 10® yr of the simulation. 



2.1 The code 

Our simulations were performed with Flash 2.2, an Adap- 
tive Mesh Refinement (AMR) hydrodynamical code devel- 
oped and made public by the ASCII Center at the University 
of Chicago ( Frvxoll ct al. 200(1) . Flash is a modular block 
structured AMR code, parallelised using the Message Pass- 
ing Interface (MPI) library. 

Flash solves the Riemann problem on a Carte- 
sian g rid using the P icccwise-Parabolic Method 
(PPM; IWoodward fc ColeUa 1 984) . It uses a crite- 
rion based on the dimensionless second derivative 

= \{Fd?F/dx^)/{dF/dxf\ of a fluid variable F to 
increase the resolution adaptively whenever T)^ > C2 and 
de-refine the grid when < ci, where ci,2 are tolerance 
parameters. When a region requires refining (X>^ > C2), 
child grids with cell size half that of the parent grid are 
placed over the offending region, and the coarse solution 
is interpolated. In our simulations, we have used density 
and temperature a s the refinement fluid variable F (see 
iFrvxell et all l2000l . for more details). We have performed 
a series of simulations with a different maximum level of 
refinement to investigate the effect of numerical resolution, 
as described in Section [2.41 

Flash interpolates the imposed analytic gravitational 
potential (see Section 12.211 to the grid cells, and computes 
the corresponding gravitational acceleration, neglecting the 
self-gravity of the gas. To interpolate the initial density field 
from our analytical solution to the Flash mesh, we impose 
an initial grid with increasing resolution from large radii to 
the centre of the cluster. This is the minimum refinement 
level allowed during the simulation. 

For the boundary conditions, we impose that the values 
of density, temperature and velocity on the boundary cells 
(guard cells) remain the ones computed at the initial time 
and satisfy the hydrostatic equilibrium conditions. With the 
cluster located at the centre of the box, the box size of 
6 X 10'^'* cm (1.9 Mpc) is large enough to ensure that cooling 
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Figure 1. The evolution of tlie total energy of each simulation is shown by plotting AE{t) = StC*) ~ £'t(0), where Eq^(t) is the sum of 
internal, kinetic and potential energy at time t and E^iO) is its initial value. The saw-tooth shape of the curves results from the discrete 
AGN events and subsequent cooling. At the mean injection rate of simulation S2.0, the energy keeps an almost constant value within 
the simulation time. We tested that this behaviour is maintained up to 5 Gyr, though we show the evolution only to 3 Gyr for clarity. 



does not affect the boundaries. This guarantees the hydro- 
static equilibrium condition is not broken, and no evident 
inflow or outflow occurs at the border of the grid. Moreover, 
the temperature difference between the last computational 
cell and the guard cell next to it, at the end of the simulation, 
is always less that 2%. Likewise, sound waves propagating 
across the volume do not generate any unphysical artifacts 
at the boundary. 

The time step is derived by the Courant condition dt = 
CAx/cs, where Ax is the dimension of a cell, Cs is the sound 
speed in that cell and C is a coefficient, usually less than one. 
The time step is uniform; thus, all the cells evolve at each 
time step. 

2.2 Initial Conditions 

The time-independent gravitatio nal potential of the dark 
matte r is th at corresponding to a iNavarro. Frenk. fc Whit3 
Jl996lll997D density profile 

Pdm/p crit 



(r/rs)(l + r/rs)2 



where pait is the critical density, rs — r2oo/c and r2oo ~ 
1.38 Mpc. We adopted a halo mass of M200 = 3 x 
10^"* MfT^ and set the concentration parameter to c 



6 iEke. Navarro fc Stei nmetz 2001; Wcchsl cr et al.l l2002t 
iDolag et al.ll2004F; ' 

W e use the prescription of IWu. Fabian, fc NulsenI 
i2000l) to set up an isothermal gas distribution in hydrostatic 
equilibrium within this dark matter potential. We adopted 
this prescriptio n because of its simplicit y. Fits to numeri- 
cal simulations iKomatsu fc Seliakll200lh . suggest that the 
temperature should decline with increasing radius. How- 
ever, our isothermal profiles do match the observed tem- 
perature structure of clusters within 0.2r2oo, where the ob- 
serv ational data show an isothermal temperature profile 
fee. lAUen. Schmidt, fc Fabianll200ll: iDe Grandi fc Molendil 
I2OOII: iPointecouteau et al.ll2004^ sometimes with a central 
decrement. A central decrement is quickly generated as the 
model cluster cools. 



We choose an initial gas temperature of T = 3.1 keV, 
consistent with the observed correlation between tempera- 
ture and M200 (RciDrich & Bohringor .20021 . The appropri- 
ate gas density is then 

PGAs/po = (l + r/r,)'"'^'-/'-=', 

where po was chosen to satisfy the relation 



,,200 /,^200 

jwgas/-'Wdm 



Here, Mp'J^ denotes dark matter mass contained within a 
sphere whose mean interior mass density is 200 times the 
critical density. We used rj = 10.25, Q — 0.3 and Q,b = 0.04 
for the matter and baryon densities in units of pcrit, and 
h = 0.7 for the Hubble parameter. 

We adopt the cooling functions of lTheuns et al.l i2002l) . 
which include cooling by H, He and metals, in the pres- 
ence of an ionising UV-background. This cooling function 
uses interpolation tables for the cooling rate due to a solar 
admixture of meta ls, obtained from Cloudy (version 94, 
iFerland et alJll998ll. We assume a met allicity of one third 
solar, and the Haardt fc Madaul ll99d) UV-backgr o und a t 
redshift z = 0, as updated by iHaardt fc Madaul i200ir) . 
The cooling routine wa s tested against the MEKAL 
(Mcwo, Lemo n, fc van den Oord 1986.,) cooling routin es in 
the X-ray analysis package XS PEC l|Kaa stra 1992) and 
against the cooling rates used bv lMcCarthv et al.. (.20041) . 

Heating by an AGN is modelled by injecting energy at 
predefined points in space and time. This allows us to simu- 
late multiple episodes of AGN activity as we now describe. 



2.3 Energy Injection 

To simulate the effect of a central AGN, bubbles of energy 
are distributed randomly inside a radius of 50 kpc from the 
centre. We do not attempt to model the processes leading to 
the formation of the bubbles, but we consider it likely that 
the bubbles will be injected at different positions because of 
the precessional motion of the jet axis of the central source. 
Energy is injected with a three-dimensional Gaussian 
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distribution proportional to e'^^ ■'^^'^ and we take the 

initial size of the bubble to be characterised by a = 10.3 
kpc. The bubble is truncated at a maximum radius of rb — 
5(7 = 51.5 kpc, which is seven cells in our standard six level 
resolution runs. The final radius (after the expansion of the 
bubble) depends on the amount of energy injected. 

Starting from t — 0, energy is injected every 10* yr 
for a time of 10^ yr in our standard simulations. We as- 
sume energy values from 1 x 10^^ to 3 x 10*''^ erg (Table 0. 
Our choice of parameters deliberately lies at the upper end 
of observed AGN duty cycles and jot powers fWillott et al. 
[1999; O wen. Eilck. fc Kassim. 2000.: Fabian ct al,..2003; Eilok. 

2.4 The Simulations 

Our main results are derived from a series of nine simulations 
(Table Q: 

1. an adiabatic test simulation (A) without cooling or 
AGN bubbles; 

2. a simulation with cooling but no AGN bubbles (SO.O); 

3. seven simulations with cooling and AGN bubbles of 
varying strength (S0.1-S3.0), but the same duty cycle. In 
these simulations only the energy of each bubble varies. The 
position and timing of the heating events are the same. 

The results of these simulations are presented in Section |21 
We also performed simulations with the following 
changes in order to test the sensitivity of our results to the 
details of the heating cycle: 

4. the bubbles were randomly distributed inside a sphere 
of 25 or 100 kpc, instead of 50 kpc; 

5. the same total number of bubbles (15) was generated 
within 1.5 Gyr, but at random time intervals; 

6. the energy was injected as bubble pairs, with each bub- 
ble containing half the energy per event. 

Finally, we also reran simulation S2.0 with increased resolu- 
tion to check for numerical convergence. The results of this 
test, described in Appendix^ led us to adopt a maximum 
refinement level of six. We also continued the simulation 
S2.0 to 5 Gyr to test the long term stability of this simu- 
lated cluster. 

In the adiabatic simulation, the gas remains almost 
static as expected for our equilibrium set-up, apart from 
small re-adjustments induced by discretisation. Energy and 
mass are conserved to ~ 0.4% and 1% respectively. With 
cooling allowed (simulation SO.O), the loss of pressure sup- 
port causes a cooling flow to be established. 

The intial conditions have a relatively high X-ray lu- 
minosity (log Lx,boi = 45.0) compared to the observed X- 
ray luminosities o f clusters of this mass and temperature 
jMarkevitchI Il998l : iMcCarthv et alJl2004 . This is a well- 
known problem that results from the low entropy gas in 
the centre of a cuspy dark matter potential. Starting with 
such a luminous cluster ensures that our simulations conser- 
vatively explore the maximum energy necessary to stabilise 
cooling. In order to be fully successful, however, our energy 
injection must not only balance the cooling, but also raise 
the entropy of the cluster gas sufficiently to reduce the total 
luminosity by a factor ~ 10. 



By default, we ran each simulation for 1.5 Gyr. With- 
out a cosmological context, a longer simulation is not well 
justifled since cluster-cluster mergers provide an important 
additional mechanism for re-organising the gas distribution 
and erasing the cooling flow structure. Nevertheless, we ran 
the S2.0 simulation for 5 Gyr to check the evolution of the 
solution over very long time scales. 



3 RESULTS 

3.1 The Effect of Energy Injection 

In Figure we show, for each simulation, the variation in 
the total energy within the simulation volume. The pulsed 
injection of the energy and subsequent cooling gives rise to 
a saw-tooth pattern. 

In simulation SO.O, which has no heating, the cooling 
rate increases with time as the cluster collapses. Including 
energy injection (simulations S0.1-S3.0) reduces the overall 
amount of energy that is radiated. This is not a trivial result 
since the energy injection events can actually promote cool- 
ing in the compressed regions around the rising bubbles. In 
practice, however, the dominant effect is the mixing of high 
and low entropy material in the large-scale gas motions in- 
duced by the bubbles. Therefore, as the injected energy is 
increased, the amount of radiated energy drops. At a mean 
energy injection rate of 6.3 x 10^^ erg s~^ (simulation S2.0), 
the radiated energy and the injected energy are balanced. 
This balance does not guarantee that the structure of the 
cluster is stable, as continuing collapse of the central region 
might be balanced by expansion of the outer part of the clus- 
ter. However, as Figure |5| shows, the whole entropy profile 
of the cluster changes little. We investigated the long term 
behaviour of this model by continuing this simulation for 
5 Gyr (Figure and found that the cluster maintains this 
behaviour even over this longer timescale. Of course, over 
such long timescales, the effects of cluster mergers can not 
be neglected. 

The emissivity-weighted temperature and entropy pro- 
files of each simulation are shown, for 50 output times, in 
Figure|21 We define the entropy as A" = kT^n^^^^ [keV cm^], 
where Te and rie are the electron temperature and density, 
respectively. This is related to the thermodynamic entropy 
by a logarithm and additive constant. If cooling is not bal- 
anced by the energy injection, the entropy profile drops dra- 
matically in the centre (top-right plot) . This is reflected as a 
marked dip in the temperature profile of the cluster (top-left 
plot). Energy injection reduces this trend, maintaining the 
relatively fiat entropy profile and isothermal temperature of 
the initial cluster. 

The cooling of material in the central cluster regions 
leads to a net inflow of gas. In Figure |3 we plot the mass ac- 
cumulation rate within the central regions of our simulated 
clusters. The rate is determined by averaging over the last 
0.3 Gyr of the simulation, and is computed within the central 
100 kpc and 50 kpc. The results for both radii are similar, as 
expected if the cluster is in a quasi-steady state. The figure 
also shows the mass accumulation rate of material with tem- 
peratures less than 1 keV. In the absence of energy injection, 
the cooling cluster deposits material in the cluster centre at 
a rate that exceeds 1000 MQyr~^, a factor ~ 10-100 higher 
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Figure 2. The emissivity-weighted temperature (left panels) and entropy (right panels) profiles of each simulated cluster are shown for 
50 simulation times. From top to bottom, simulations SO.O (without any AGN bubble), SO. 3, Sl.O and S2.0 are shown. The effect of 
energy injection is to maintain the initial, relatively shallow entropy profile and isothermal temperature in the cluster core. 



than the observational limits JCrawford et alJ 119991 : lEded 
I2OOII: ISalome fc Combe3l2003l) . Similarly, this model has a 
large mass flux of material cooling below 1 keV (1/3 of the 
ambient temperature). 

Energy input at an average rate of 4.7 x 10** erg s~^ 
(SI. 5) reduces the deposition rate below 30 MQyr~^ at 
the end of the simu l ation, as required by observations 
fPctcrson ct al." '200 1'; 'Ta mura et"al] 1200 ll: iKaastra et"ai] 
I2OOI : Peterson et al. .2003,) . Even with a lower energy in- 
jection rate of 3.2 x 10** erg (Sl.O), the rate at which 
material cools below 1 keV is still consistent with observed 
X-ray spectra. We conclude that energy injection through 
the creation of buoyant bubbles can successfully prevent 
catastrophic cooling. 

3.2 Two-dimensional morphologies 

In Figure |1] we show the time evolution of various quan- 
tities in two-dimensional sections of simulation S2.0. The 



temperature and entropy distributions are shown for a slice 
in the x-y plane. They show that cold, low entropy gas is 
pushed to large radii by the bubble's buoyancy and mixed 
with gas at those radii. Thus, the cooling rate is reduced 
not by direct heating of the cooling gas, but by convective 
transport of this material to regions of lower pressure. Sim- 
ilar results have been seen in simu l ations of single bubbl es 
(e.g. iQuilis. Bower, fc Baloghll200lt IChurazov etalll2002h . 

The temperature distribution also reveals the presence 
of sound waves propagating through the ICM. The sound 
waves are almost concentric and regular, a consequence of 
the periodic energy injection events near the cluster centre. 

In order to illustrate what might be observable with 
an idealised X-ray satellite, we also calculated the projec- 
tion of the quantity p^T^^^, which is approximately propor- 
tional to the bolometric gas emissivity. Despite the large 
amounts of energy being dumped into the central regions, 
the large-scale gas distribution appears smooth and undis- 
turbed in projection. To reveal the presence of perturbations 
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Figure 4. The time evolution (from left to right) of simulation S2.0. From top to bottom, the quantities shown arc the temperature [K] 



and entropy [ergs g 



-2/3 



cm^] on the x-y plane crossing the centre. In the bottom row we show the approximate bolometric emissivity 



p2 jil/2 projected through the simulation volume. The temperature distribution reveals the presence of sound waves propagating through 
the ICM. The sound waves are almost concentric and regular, a consequence of the periodic energy injection events near the cluster 
centre. 



on this smooth profile, we created an unsharp masked image 
through the transformation (imgg — imgg)/imgj,, where imgg 
and imgj, are the original image and its smoothed version, 
respectively. The initial image was taken from a high resolu- 
tion simulation (~ 2 kpc) and the final unsharped image is 
shown in Figure|3 The contrast in this image is high enough 
to show the sound wav es distinctly as ripples. As suggested 
bv lFabian et alJ ll2003l) . these ripples in X-ray images of ob- 
served clusters could be the fingerprint of multiple outbursts 



and could be used to measure their periodicity. Furthermore, 
if the cluster gas is sufficiently viscous, these sound waves 
mig ht help to offset the cooling fiow by directly heating the 
gas (iFabian et al]l2003l : iRuszkowski. Brfiggen. fc Beeelmanl 
200^. By comparing the properties of successive ripples it 
may be possible to measure the viscosity of the ICM. Al- 
though the ripples in our simulation extend to very large 
radii, a result of the small viscosity in the simulations, it 
is unlikely that realistic observations would have sufficient 
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Figure 3. Top panel: The mass deposition rate into a sphere 
of radius 50 kpc (dashed line) and 100 kpc (dash dotted line), 
averaged over the time interval 1.2 to 1.5 Gyr, as a function of the 
energy injected per bubble. For bubbles of energy above 10®" erg 
the deposition rate is < 100 M0yr^^. Bottom panel; The mean 
growth rate of the total mass at temperature < 1 keV for each 
simulation. The amount of gas cooling below 1 keV drops down to 
almost zero when the injected energy per bubble is greater than 
10'^° erg. 



signal-to-noise ratios to see these weak features at such dis- 
tances. We will investigate this further in a future paper, 
where we synthesise artificial X-ray images, with the appro- 
priate transmission functions, from our simulations. 



3.3 Sensitivity to duty cycle parameters 

Simulations S0.1-S3.0 were based on a regular injection 
of energy. This is somewhat artificial, so we investigated 
whether randomly-timed energy injection events would have 
the same effect. We also investigated what happens if the 
energy injection occurs in less frequent bursts of greater in- 
tensity, or if bubbles are distributed inside a sphere of larger 
(100 kpc instead of 50 kpc) or smaller (25 kpc) radius. All 
our tests were performed choosing the parameters of simu- 
lation S2.0. 

The resulting energy evolution of each simulation is il- 
lustrated in Figure |n] First, we double the time between 
bursts, while keeping the total average energy rate constant; 
thus, each bubble has twice as much energy injection as in 



S2.0. In this case, the variation in total energy during each 
cycle is much larger than before. However, the long term 
evolution is very similar. Randomising the injection time 
does not have a large effect on the average energy evolution 
over the duration of the simulation either. Where a particu- 
larly long time interval elapses between bubbles, the energy 
drops further as cooling progresses; however this is offset 
by later events where bubbles appear in quick succession. 
Our results are also insensitive to whether the bubbles are 
injected singly, or in pairs with the total energy divided be- 
tween the two bubbles, possibly a better model for the AGN 
activity. The main factor that does make a significant dif- 
ference is the location of the bubbles. Placing the bubbles 
initially within a larger radius (100 kpc) is much less effec- 
tive at disrupting the cooling flow, since these bubbles tend 
to rise buoyantly without disturbing the central, cooling gas. 
By the end of this simulation, the total cluster energy has 
dropped below that of simulation SI. 5 and is falling rapidly. 
Conversely, placing all the bubbles within 25 kpc is much 
more effective at preventing cooling, and the total energy 
of the system actually increases with time. To summarise, 
we see that the energy balance is most sensitive to the total 
amount of injected energy, and to the initial location of the 
bubbles. 



4 DISCUSSION 

4.1 Energy Requirements 

Most central galaxie s in clusters are radio sources 
iLedlow fc Qwenll99g ) ; we can therefore expect that events 
like those we have simulated will play an important role 
in cluster evolution. We have shown that the cluster cooling 
rate can be reduced so that it agrees with observations if the 
AGN activity is sufficiently energetic. We have modelled the 
energy injection events as short outbursts separated by rel- 
atively long quiescent periods in order to make an extreme 
test of the long term effect of the energy injection mecha- 
nism. To stabilise the long term evolution of the model clus- 
ter in this way, we require an energy of ~ 2 x 10®" erg for 
each of the injection events, if the bubbles are distributed 
within a ~ 50 kpc radius. We can associate these events 
with outbursts of powerful AGN activity with a duration of 
~ 10^ yr and a power of ~ 6 x 10*^ erg s~^. 

The observed radio power of central galaxies varies by 
factors of ~ 100 for clusters with similar core X-ray luminosi- 
ties. The most powerful radio sources can have monochro- 
matic radio power of u ■ Sv ~ 10^^'^ erg s~^ (where Si, is the 
radio luminosity density a,t i/ — 1.4 GHz); around 10% of 
clusters host such powerful sources. Eilek's (2003) analysis 
of central cluster radio source morphologies and ages sug- 
gests that the ratio of the total jet power (Pjct) to u ■ Si, is 
more than 10^. Under minimal assumptions, this can be con- 
firmed in the case of M8 7, where the most detaile d radio ob- 
servations are available JOwen. Eilek. fc Kassirnf2 000 ) ; such 
pow erful jets are also inferred for classical FR H sources 
iWiUott et al.iri99a.) . This power is close to the Eddington 



limit for accretion by a 10 Mq black hole. Thus, the total 
jet power of the most powerful radio galaxies is compara- 
ble to the energy required to counter-balance cooling in our 
simulated cluster. The frequency with which these sources 
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Figure 5. An unsharp masked projection of tiie approximated 
X-ray emissivity, p-^T^/-^, for the central (0.9 Mpc diameter) part 
of tlie cluster simulated at a maximum resolution of ~ 2 kpc. 
The "granular" structure of the image is due to the resolution of 
the simulation grid. The ripples from successive energy injection 
events can be clearly seen. 



occur in clusters also seems consistent with the duty cycle 
that we have assumed in our simulations. 

It is worthwhile to emphasise that the energy injection 
can vary significantly (by a ~ 50%) while still reducing the 
mass deposition rate to an acceptable level after 1.5 Gyr. 
These variations in the injected energy lead to an overall 
cooling or heating of the cluster: however subsequent merg- 
ers with massive substructures will disrupt the evolution 
presented in our simulations (which assume a static grav- 
itational potential). The shocks and turbulence generated 
in the mergers will lead to a re-distribution of energy and 
entropy that may erase the differences that have built up 
during the quiescent phase that we simulate llGomez et alJ 
|2002; Burns ct al. 2003). Our model does not contain any 
mechanism to regulate the energy injection - this would re- 
quire detailed understanding of the physics powering the 
ce ntral source from the su rrounding gas reservoir. However, 
as iBinnev fc TaboJ (|l993) (see also iSinne^liooi) have ar- 
gued, it seems natural that the cooling rate and the energy 
injected by the central energy source should be linked. In 
this balance between heating and cooling is nat- 

urally expected when these quantities are integrated over 
sufficiently long timescales: if insufficient energy is injected 
into the ICM, the net cooling leads to a high mass deposition 
rate, and thus an increase in the fueling of the central en- 
gine. Nevertheless, the delays inherent in this feedback loop 
lead to a short-term imbalance, and thus to the out bursts 
that we model here. 



4.2 Global X-ray properties 

Although the injected energy is able to counterbalance the 
cooling flow in the cluster, the effect on the cluster's overall 
structure is modest, as we showed in Figure |21 In the ab- 
sence of any heating, cooling causes the central entropy and 
temperature to drop precipitously. In simulation Sl.O, the 
AGN energy is sufficient to reduce this drop to only ~ 50 
per cent. Thus, by preventing runaway cooling, our energy 
injection mechanism reduces the change in the initial en- 
tropy profile of the cluster. The result is that the integrated 
luminosity and average temperature of the cluster evolve lit- 
tle. Doubling the amount of energy injected (S2.0) results 
in a small increase in the central temperature (~ 30%), and 
a 20% reduction in luminosity. 

We simulated an initially isothermal cluster, with 
temperature observationall y consistent with its mass 
jReiprich fc Bo hringcr '2002fl. in order to see how it would 
evlove in the luminosity-temperature plane. The experi- 
ment raises two issues that are not fully resolved by our 
model: (1) The cluster is too luminous, for its emission 
weigh ted temperature, compared with observed clusters 
(e.g., 'Markcvitch 1998^ , even w hen uncorrected for cool- 
ing flows (McCarthv et alj|2004l) . (2) The temperature pro- 
files of the simulated clusters do not show the temperature 
decrement seen in the central regions of real clusters (e.g. 
^Ucn. Schmidt, & Fabian 2001). 

In order to lower the luminosity of the simulated clus- 
ter enough to agree with observations, we can either ad- 
just the gravitational potential (by assuming a lower halo 
concentration) or increase the entropy of the central gas. 
In order to get good agreement with the average lumi- 
nosities of clusters of this temperature, we would need to 
adopt c ~ 3, well below the concentration s inferred from 
dark matter simulat ions {Eke. Navarro, fc* S tcinmctz 200^; 
IWechsler et al.ll2002l) . On the other han d, the central cn- 
tropy could be raised to > 160 keV c m~'^ jBabul et al. 200^ : 
IVoit. Brvan. Baloeh. fc 13oweil20o3l , which is a factor of five 
greater than the initial entropy at 20 kpc. Even if enough 
energy were injected to achieve this (for example by run- 
ning simulation S2.0 for a much longer time), this would 
make the core temperature even hotter, exacerbating the 
disagreement with observed temperature profiles. Thus, we 
conclude that the initial cluster gas distribution must have 
been different from the isothermal model we have assumed. 
Since the high cooling rate of the unperturbed cluster is a 
direct consequence of its low central entropy (leading to high 
X-ray luminosity), our energy injection mechanism has re- 
moved the symptoms of the overcooling problem, but not 
addressed its cause. 

Global preheating of the intragalactic medium 
prior to cluster formation should lead to clus- 
ter gas distributions with integrated luminosities 
and temper atures t hat are consistent with ob- 
serva tions iKaisei] |l99l|: [Navarro^ Frenk. fc White ! 

'i895j ^ 

Balogh. Babul, fc PattonI 119991 : iMuanwong et al.1 
Borgani et aL, . 2_00 j) . An alternative (but related) possi- 
bility is that heating events similar to the ones we have 
simulated here occur in the cluster progenitors and/or 
surrounding filamentary structure, at earlier times. Raising 
the entropy in these lower mass structures may result in 
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Figure 6. In each panel we show the evolution of the total energy in simulations S2.0 and SI. 5 (dotted lines), as in Figure □] The other 
lines represent models with the same total energy as simulation S2.0, but with different parameters for the energy injection mechanism. 
In panel (a) we show a simulation where the duty cycle is regular, but twice as long as in S2.0. In panel (b) we show the effect of 
randomising the time interval between injection events. In panel (c) we show the effect of injecting the energy within pairs of bubbles. 
Finally, in panel (d) we show, as the dashed and continuous curves, the effect of placing the bubbles randomly within a 25 kpc or 100 
kpc area, respectively. 



a higher entropy for the final cluster, as the lower density 
results in an increase in the accret ion shock strength 
JVoit et al.ll2003l : IVoit fc Ponmanll2003l) . In either scenario, 
subsequent cooling would naturally lead to a temperature 
decrement i n qualitative agre e ment with observations. Cal- 
culations by iMcCarthv et al.l (|200^ suggest that although 
the cluster would then increase in luminosity, it would 
decrease in emission weighted temperature so as to remain 
in acceptable agreement with the observed L-T relation. 

4.3 Entropy transfer and generation 

In Figure we have shown that an energy injection rate 
of 6.3 X 10** erg s~^ is sufficient for the total energy of the 
cluster to remain constant and for the density and temper- 
ature distribution of the ICM to reach a quasi-equilibrium 
profile in which the radiative heat losses are balanced by the 
energy deposited in the ICM by the hot bubbles. The slow 
evolution of the density profile implies that the energy we 
inject is somehow shared with the cooling ICM and does not 
remain trapped in the hot material. In particular, the energy 
must be efficiently shared with the dense material near the 



centre of the cluster, where the cooling time is short, leading 
to an overall increase in entropy in these regions. 

There are three mechanisms by which the energy can 
be shared: 

(i) PdV work - as energy is injected into the bubble, it ex- 
pands doing work on the surrounding material. We calculate 
the PdV work by calculating the volume of the bubble at 
a series of timesteps. 40% of the injected energy is used in 
this way. The process is largely adiabatic, however, and has 
little efi'ect on the entropy distribution of the surrounding 
ICM. 

(ii) shocks - If the motion of the the bubble were super- 
sonic, shocks would dissipated energy in the surrounding 
ICM leading to an overall increase in the entropy of the 
system. In practice, however, the bubble does not reach su- 
personic speeds: as the speed of the bubble increases the 
energy dissipation in the turbulent wake grows. In our sim- 
ulations the maximum Mach number in the flow is 0.8 thus 
we do not expect shocks to be a major source of entropy in 
the surrounding ICM. 

(iii) up-Iift and turbulence - The rising bubble generates 
a complex flow pattern, with material from the central re- 
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gions of the cluster being up-lifted in the wake behind the 
bubble. Vortices in this flow dissipate energy and mix the 
bubble material, low entropy material from the center and 
the surounding ICM. These irreversible processes generate 
an increase in the entropy of the cluster. 

Even though we assume that the viscosity of the ICM 
is small, it is certainly not negligible, and viscous dissi- 
pation on molecular scales can still be a major source of 
heating in turbulent regions of the flow. Clearly our simu- 
lations are far from resolving the relevant molecular diffu- 
sion scales. However, in turbulent flows, the non-linear cou- 
plings ca use a cascade of e nergy from large to small scale 
motions jKolmogorovlll94lfl . The viscosity determines the 
smallest scales in the flow, while the rate at which energy 
is dissipated is determined by the largest scale eddies: these 
control how quickly energy is fed into the turbulent cascade 
(Tcnnokcs & Lumlcv 1972). In simulations, the development 
of the turbulent cascade is limited by the numerical mixing 
of fluid elements within a grid cell, rather then molecular 
viscosity. However, the large scale properties of the flow are 
not expected to be dependent on the smallest scales that 
are resolved. We check this by comparing the evolution of 
a single bubble in two simulations: namely the fiducial res- 
olution and at almost an order of magnitude greater linear 
resolution. Plotting the two simulations to the same base 
resolution we see that the gross structure of the flow is the 
same. Thus, because of the weak dependence of the large 
scale structure of the flow on the numerical resolution, it 
is plausible that the energy dissipation and mixing rates in 
the regions around the rising bubble will converge over the 
range of resolution that we can test in Appendix 1X1 

In addition to being a source of dissipative energy, the 
turbulence in the flow provides a means of increasing the 
entropy of the surrounding ICM through mixing. The mix- 
ing may occur between the hot material within the bubble 
and the surroundings, or between the surrounding ICM and 
the low entropy material drawn out of the cluster centre by 
the up-draft created by the rising bubble. In both cases, the 
mixing is irreversible and leads to an overall increase in the 
entropy of the system. In the first case, the mixing transfers 
energy from the bubble to the surrounding ICM (in addition 
to the PdV work discussed above) . The second case does not 
directly transfer energy from the bubble, but by driving an 
entropy increase in the lowest entropy material (the mate- 
rial with the shortest cooling time) it prevents this material 
cooling out of the flow. 



5 CONCLUSIONS 

We have presented 3-dimensional gasdynamical simulations 
using the Flash adaptive-mesh refinement code of a cool- 
ing flow in an isothermal. X-ray luminous cluster of galax- 
ies, with periodic injections of thermal energy. The initial 
cluster has mass of 3 x 10^* Mq and a temperature of 
3.1 ke V, consistent with the observ ed mass-temperature re- 
lation l|R,eiprich fc BohringeilEoO^ . The injected energy is 
manifested as hot bubbles of buoyant gas that rise out of 
the cluster core, convectively mixing the cooling material. 
We treat these injection events as sporadic outbursts with a 
typical duty cycle of 100 Myr. The parameters of these in- 
jection events match the observed luminosity and frequency 



with which the most powerful radio galaxies are seen in clus- 
ters. 

In the absence of any energy injection, mass rapidly 
flows to the cluster centre, at a rate that exceeds observa- 
tional limits by at least an order of magnitude. Based on 
simulations with a variety of parameters related to the en- 
ergy injection rate, run for between 1.5 Gyr and 5 Gyr, we 
draw the following conclusions: 

1. For a time averaged energy injection rate of 6 x 
10*'' erg s~^, the mass inflow rate is less than 30 MQyr~^, 
compatible with available observational limits. With this 
amount of heat input, the total energy of the cluster remains 
approximately constant over 5 Gyr. 

2. The evolution of the total cluster energy depends pri- 
marily on the total amount of energy injected and the spatial 
distribution of bubbles, but is only weakly sensitive to the 
duty cycle of heating events or to whether the bubbles are 
produced singly or in pairs. 

3. The bubble activity generates concentric sound waves 
that are clearly evident in unsharp-masked projections of 
cluster emissivity. 

4. When the injected energy just balances cooling, the en- 
tropy and temperature profile of the cluster remain approx- 
imately unchanged from their initial configurations. 

In summary, periodic energetic events of the kind we 
have simulated can reduce the mass flow rate and accumu- 
lation of cold gas in massive clusters to within observational 
limits. However, this mechanism operating on a fully formed 
cluster does not result in a flnal luminosity consistent with 
observations. It is likely that the structure of the progeni- 
tors from which the cluster formed was affected by heating 
events prior to the assembly of the final cluster. 
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Figure Al. As Fig. the energy evolution is shown for the 
simulations Sl.O, SI. 5 and S2.0 run with 6 levels of refinment 
(black, solid lines). We compare these with simulation S2.0 run 
at different resolutions (7, 8 and 9 refinement levels). 



APPENDIX A: RESOLUTION CONVERGENCE 
TESTS 

We tested the convergence of the code by running simula- 
tion S2.0 at increased resolution. In Figure IXTl we compare 
the total energy evolution for simulation Sl.O (6 levels of 
refinement, resolution of 7.6 kpc) and S2.0 (6, 7 (3.8 kpc), 
8 (1.9 kpc) and 9 (0.9 kpc) levels of refinement). The codes 
were run for different lengths of time because of the lim- 
itation of the computational time and memory available. 
The differences between the four versions of the S2.0 run 
are small in comparison to the relative evolution of the SI. 5 
and S2.0 simulations. The figure also shows that the radi- 
ated energy does not increase systematically with resolution. 
Comparing the total radiated energy at t = 0.45 Gyr, the 
four versions of the S2.0 simulation differ by less than 3%. 
At the last common output time, the 6 and 7 level sim- 
ulations differ by 4%. We also compared the density and 
entropy profiles of the simulations at this output time, and 
found a similar level of agreement. We therefore chose to run 
the code with a maximum refinement level of 6, this being 
a good compromise between the speed of the code and the 
desired accuracy of the relevant quantities. 

Simulation SO.O was also used to check the convergence 
of energy conservation as the maximum number of levels of 
refinement changes. With just 4 refinement levels, the total 
energy has an error of ~ 1%, while with 6 levels the error 
is less than 0.1%. We decided to run the simulations with a 
maximum number of 6 levels. The equivalent resolution of a 
fixed-grid Eulerian code is 256'^. 



